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SUMMARY 


An iterative numerical method has been developed for the 
calculation of steady, three-dimensional, viscous, compressible flow 
fields in centrifugal compressor impellers. The computer code, which 
embodies the method, solves the steady three-dimensional, compressible 
Navier-stokes equations in rotating, curvilinear coordinates. The 
solution takes place on blade-to-blade surfaces of revolution which 
move from the hub to the shroud during each iteration. 

Numerical calculations were made for two centrifugal impl- 
lers, one with radial blades and the other with backswept blades. The 
radial impeller operated in a laminar Reynolds number regime. The 
backswept impeller problem was used to check out the turbulence model 
incorporated in the code. A large vortex was calculated on the suction 
blade surface of the radial impeller in the region of the discharge; 
such a vortex is qualitatively in agreement with observations. The 
backswept rotor calculation did not indicate an impeller separation. 

No conclusions can be drawn with regard to the effectiveness of backsweep 
in reducing or eliminating flow separation, because the radial impeller 
was calculated for laminar flow at a very low Reynolds number (5000) , 
whereas the backswept impeller problem was calculated for turbulent flow 
conditions with the turbulence model operational. Contour plots are 
presented to show the calculated static pressure in the blade-to-blade 
channel. Relative velocity vector plots on various blade-to-blade surfaces 
show significant differences with inviscid potential flow solutions in 
common industry usage for centrifugal compressor design. It is concluded 
the viscous Navier-Stokes solution for flow fields in centrifugal compres- 
sors represents a significant advanceiuent in the ability to analyze these 
complex types of turbomachinery. 


1.0 


INTRODUCTION 


The principal objective of this research effort is to develop 
a computer program to calculate the three-dimensional, viscous, 
compressible flow field in blading passages of turbomachinery . At present 
the computer program is being applied to the rotating impeller of centrifugal 
compressors. 

This submittal is a report on the work which has been completed 
in Phase I of a two phase program of research and development. Three 
main tasks have been completed in Phase I. They are as follows; 

1. An impeller computer code was developed and debugged. 

2. A radial centrifugal impeller problem was solved. 

3. A backswept centrifugal impeller problem was solved. 

Phase II of this research effort is comprised of two additional principal 
tasks . 

4. To speed-up the computer code of Phase I by a factor 
between 3 and 5. 

5. To revise the computer code of Phase I, which calculates 
the flow field on blade-to-blade surfaces, to calculate 
the flow field on cross-sectional surfaces. 

The importance of the cross-sectional calculation is discussed in Section 

4.0. 
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2.0 SYMBOLS 


Cy 

E 

H 

hx 

£ 

io. 

^2 

-i-3 

2 

J 

K 

k 

Ks 

M 

m 

n 

P 

r 





t 

u 

u 


Uoo 

V 

W 

w'=w-u 

X 


Heat Capacity at Constant Pressure 

Heat Capacity at Constant Volume 

Specific Internal Energy 

Thermodynamic Heat Function or Enthalpy 

Metrics of Transformation 

Metrics of Tramsformation 

Metrics of Transformation 

Unit Vector of Curvilinear Coordinate x 

Unit Vector of Rotating Cartesian Coordinate of 

Unit Vector of Rotating Cartesian Coordinate of X2 

Unit Vector of Rotating Cartesian Coordinate of X3 

Unit Vector of Curvilinear Coordinate y 

Index Specifying Streamlike-lines on blade-to-blada Surface 
Index Specifying Potential-like-lines on Blade-to-blade surface 
Unit Vector of Curvi liner Coordinate z 

Von XartTian's Constant 

Momentum 

Mass 

Time Index for Finite Difference equation 
Pressure 

Maxium Radius of the Impeller (at the exit) 

Radial Coordinate which together with X3form a Cylindrical 
Coordinate System 

Grid Velocity Component along x Direction 
Grid Velocity Component along y Direction 
Total Laminar Stress Tensor 
Time 

Particle Velocity Component along x Direc**ion 

Particle Velocity Vector 

Speed of March along z Direction 

Particle Velocity Component along y Direction 

Particle Velocity component along z Direction 

QQ Velocity along z on a Galilean Frame which ireves with a 

Constant Speed Uqo along z with respect to the laboratory frame 

Curvilinear Coordinate along Azimutlial Direction 
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Coordinate Axes of Rotating Cartesian Coordinate which Rotate 
about Axial Axis X3 with Speed w 

X2 Coordinate Axes of Rotating Cartesian Coordinate which Rotate 

about Axial Axis X3 v ith Speed w 

X^ Axial Coordinate 

Y Curvilinear Coordinate along Streamwise Direction (from inlet 

to discharge) 

Z Curvilinear Coordinate or. Marching Direction (from h\i> to shroud) 


Y 

S 

€ 


A 

A 

A 

e 


Symbols in Greek Letters 

Heat Capacity Ratio Cp/'C^ 

Boundary Layer Thic)cness 
Incompressible Displacement Thic)tness 
Eddy Viscosity 

Molecular Viscosity Coefficient 
Kinematic Viscosity Coefficient 
Rotation Velocity of Impeller 
Total Stress Tensor 
Reynold Stress Tensor 

Pressure Blade Surface Meridional Angle 

Local Flow Angle Between Pressure Blade Surface and Meridional plane 
Density 

Shearing Stress at Wall 

Viscosity Coefficient for the DeViatoril Strain = 

Azimuthal Coordinate Angle, together wit)i r and X3 form 
cylndrical coordinate system 
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3.0 


BACKGROUND 


In recent years, considerable effort has been spent in solving 
the time -dependent, compressible, Navier-Stokes equations for systems with 
plane two-dimensional and/or axial symmetry 2, 3, 5,6,7. single 
numericad method was used to solve these two-dimensional and/or axisym- 
metric problems. The numerical method* is an explicit time marching scheme 
in two spatial dimensions. Details of the method are presented in References 
4 cind 7 . 

In 1969, under sponsorship of NASA Ames Research Center, a research 
effort was initiated to apply the above time-dependent, two-dimensional 
method to solve steady flow problems in three spatial dimensions. The basic 

Q 

idea was based on the Equivalence Principle °, which states that for 
slender bodies at hypersonic speeds the three-dimensional steady equations 
of motion for inviscid flow reduce identically to unsteady equations in two 
dimensions . 

This principle was extended in an ad hoc manner to a viscous 
flow through a model which permits viscous cross flow together with inviscid 
axial flew. Figure 1 shows an ogive-cylinder body at angle-of-attack 
with respect to the freestream flow direction; leeward vortices are also 
indicated in the figure. The axial coordinate z was made proportional to 
a time-like-variable, t, according to the relation 

z = Uoo-t (1) 

where Uco is the freestream speed. The two-dimensional Navier-Stokes 
equations were solved in cross-sectional planes normal to the system's 
axis. The cross section planes wore moved at freestream speed from the 
leading edge to tlie trailing base of the body. The time-dependent flow 
field at each cross-sectional plane corresponded to the steady flow at Uie 
z coordinate given by Equation (1). Since the leeward vortices of Fi<iure 1 
have axes which are almost parallel to the z axis, tiie cross-sectional 


* The numerical method was originally developed by Trulio. 


- 4 - 




planes of calculation contain these vortices. Although all axial effects 
were neglected, this numerical procedure did calculate leeward vortices 
amd produced other flow field results which were generally in accord with 
experiment 

Subsequent to the above research effort, a new study was 
launched, under NASA Ames Research Center sponsorship, whereby axial effects 
were incorporated in the nun;erical procedure and a numerical solution to the 
full Navier-Strokes equations wore generated by iteration. A one-dimensional, 
time-dependent radial computer code was used to solve for the steady, viscous, 
compressible, supersonic flow field about a cone-cylinder-flare body at zero 
eingle of attack. In three iterations the calculated boundary layer and shock 
wave structure converged. After five iterations a recirculation region 
formed at the cylinder-flare junction. Although the computed recirculation 
region was much smaller than measured, cemparisions of calculated shock 
structure and boundary layer results with experimental data and boundary 
layer theory predictions were satisfactory*. 

The above axisymnietric computer code was then revised to solve 
for two-dimensional time-independent flow fields. Trulio cind Yeung solved 
for the supersonic viscous, compressible flow fieid in a ramp-compression 
corner**. As in the case of the cor.c-cylinder-f lare junction, the iterative 
method did not accurately predict the recirculation region formed at the shock- 
wage-boundary layer interaction. 

The iteration procedure employed for the axisymetric and two- 
dimensional flows, described above, was completely reformulated to account 
for boundary layer separation and the subsequent evolution of a recirculation 
regim. The re-formulated iteration procedure was then appliec, to the impeller 
of a centrifugal compressor in tl;is research effort. Formulatitn of the 
iteration procedure for the impeller problem is discussed in tlie next section. 


*Walitt, L. , "Computation of Steady Axisymmotric Flow Using a One-Dimensional 
Time Dependent Method, "Applied Theory Report ATR-74-16-1, August 1974, to be 
published as a N/eSA Contractor Report. 

**Truilo, J.G., and Yeung, H.W., "Iterative Solution of the liquations tor 
Steady Viscous Compressible Flow Based on Similitude," Aerospace Rescarcli 
Laboratory, Report No. ARL74-0138, 1974. 
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4.0 FORMULATION OF IMPELLER PROBLEM 

A set of finite difference analogs of tlie full three-dimensional, 
compressible, Navier-Stokes equations has been developed and programmed. 

In addition to three-dimensionality and compressibility, the following 
effects are included: 

1. Centrifugal Force 

2. Coriolis Force 

3. Transition and Turbulence 

4. Arbitreiry Impeller Geometry 

5. Impeller Tip Clearance 

A solution to these finite difference equations is obtained in the following 
manner. For the radial impeller, eui inviscid flow field was generated by 
the method of Reference 10. For the backswept impeller, an inviscid flow 
field was generated by the method of Reference 39. Stc.rting from the known 
inviscid solution, the viscous effects are calculated through iteration. 
Certain terms of the finite difference equations (FDE) are evaluated from 
the inviscid solution and other terms are evaluated directly. Terms evaluated 
from the inviscid solution are designated "elliptic source terms", while 
those evaluated directly are designated, "parabolic terms". 

The distribution of the elliptic source terms and parabolic terms 
in the FDE depends on the mode of marching. At present two modes of marching 
are contemplated. 

1. The FDE are solved on blade-to-blade surfaces which move 
from the hub to the shroud. 

2. The FDE are solved on cross-sectional surfaces which move 
from the inducer to discharge. 

Each method of marching results in its own set of elliptic source terms and 
parabolic terms. 

For illustrative purjxjses we start with a schematic of a typical 
impeller for a centrifugal compressor shown in Figure 2. In the Llade-to- 
bladc mode of marching, the computation takes place on a blade-to-blade 
surface, which extends from inducer to tlie discharge, and moves from Uie 
hub to the shroud during an iteration. n-ie darkened surface of Figure 2 is 
the hvib blade-to-blade surface. The blade-to-blade method of marching is 
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illustrated in the blade passage schematic shown in figure 3. The x, y, and 
z coordinates of Figure 3 represent a loft handed orthogonal, curvilinear 
coordinate system. The z - direction is proportional to the timc-like- 
variablc, t, with the calculation taking place in the (x,y) Llade-to-blade 
surfaces. The (x,y) blade- to-blade surfaces move from the hub to the 
shroud of the impel le-. This mode of maurching accounts for two very important 
fluid mechanical effects that occur in impellers. 

1. Upstreaun influence effects - 

The flow is subsonic relative to the moving blidcs; hence, downstream condition 
influence upstream conditions. Since, each blade-^o-blade E .rfacc extends from 
inducer to discharge, the downstream flow cox. influence the upstream flow as 
the blade-to-blade surface moves from the huo to shroud. 

2. Blade boundary layer separation - 

Separations, which occur on the blade surfaces, produce vortices whose axes 
are normal to the blade-to-blado surfaces. Thus, the vortices themselves 
are contained in the blade- to-blade surfr.ee and are easily calculable. 

The cross-sectional mode of marching is analogous to the body-at- 
axigle-of-attack problem discussed in Section 3.0. We march down the channel, 
from the inducer to discharge, in cross-sectional surfaces normial to the hub 
surface. A schematic of the blade passage with the c.rosE-sectional surfaces 
indicated is presented in Figure 4. The z-coordinate , which varies 
with time, is now normal to the (x,y) cross-rectional surface of Figure 4. 

The (x,y) cross sectional surfaces move from the inducer to the discharge 
of the impeller. This mode of marching accounts for three additional fluid 
mechanical effects that occur in impellers. 

4. Channel corner vortices 

At the junctions of the bl.,des and the imb, vortices usually form, who.se axe.s 
are generally normal to the cross-sectior.il surface;.; hence, the corner 
vortices would be contained in these surfaces and are easily calculai.le . 

5. Shroud effects - 

Relative to the blades , tlif; shroud imposes a moving boui.dar>' conditioii. T!ie 
effects of this moving boundary condition may induce sej^aration in tlu- 
neighborhood of the shroud. This separation is calculalilc in ci o;;s-sect ional 
surfaces since each surface contains the shroud vortices. 


-'J- 



















6. Blade tip clearance effects - 

Since the shroud and blade tip are contained in each cross-sectional plane, 
spillage in the tip clearance region is calculcible in this mode of marching 

To properly solve for am impeller flow field, an iteration procedure 
with both modes of marching is required. The procedure is as follows. 

Starting from an inviscid solution as the "zeroth" iterate, we determine 
the first viscous iterate by maurching in blade-to-blade surfaces which move 
from the hub to the shroud. Based on the first iterate we determine a second 
viscous iterate by marching in cross-sectional plaines which move from the 
inducer to the discharge. In this way the six principal impeller fluid- 
mechanical effects, described abo'.e, can be accounted for. The second 
iterate will be a complete solution to the three-dimensional, compressible, 
Navier-Stokes equations for flow in a centrifugal impeller. 

The blade-to-blade mode of marching has been developed in Phase I 
of this research effort. The blade-to-blade results generated are the 
sxabject of this report. 


5.0 


DFRIVATION OF THE INTEGRAL CONTINUITY EQUATION 


In this section the integral continuity equation solved on the 
(x,y) blade-to-blade surface is derived. This derivation is presented 
to illustrate the actual elliptic source terms and the pareibolic terms 
of the equations of motion. The equations of motion in rotating, orthogonal, 
curvilinear, Eulerian coordinates x, y, and z are presented in Appendix A. 

The steady three-dimensional equations of Appendix A, in 
Eulerian coordinates x, y, z, are transformed to (x,y,t) space according 
to the following relations: 


z = Uoo t , 


A. = A. 

dz Uoo ' 


U 


oo 


+ w 


( 2 ) 


where t is a time-like-variable, U is the velocity of the blade-to-blade 
surface, w is the velocity component in the z-direction, and w’ is the 
velocity in the z-direction. Equations (2) represent a mathematically conven- 
ient trcinsformation and lead to a compact set of integral relations; however, 
they are somewhat non-physical in that the variable t may no longer be time- 
like, having the units z/U 

The conseravation of mass for steady motion relative to the rotating 
curvilinear coordinates (x,y,z) is (Appendix A' as follows: 

AS) + ^ V” AK) = 0 (3) 


where is the density, u the x-velocity component, v che y-velocity component, 

and (h ,h ,h^) are the trcinsformation metrics. 

X y t 

According to Equations (2) the continuity equation becomes: 



The left-hand-side of Equation (4) closely resembles the continuity equation 
for unsteady flow in the (x,y) plane. The transformation metrics h^, h^, 
on the left-hand-side account for the fact that the flow is not planar but 


occurs on a curvp^i surface. The term on the right-hand-side of Equation (4) 
represents a source term which accounts for the variation of axial velocity 
w from the consteint reference velocity. This term must be considered known 
in the iteration process and is evaluated from the previous iterate in 
each successive iteration. 

Equation (4) is formulated in the rotating, Eulerian coordinates 
x,y,t, however, the calculational process takes place on the (x,y) blade- 
to-blade surface (Figure 3) which distorts with time t according to the shape 
of the blade surfaces. Hence, we are really interested in the continuity 
equation in a generalized coordinate system ^ / where 

t = r 

X = f (]p , 7 , t ) (5) 

y = g if 

The transformed continuity relation in ^ , 7 » ^ space, which is derived in 
Appendix B, is presented below. 


A 


A 



C 


( 6 ) 


where 




(7) 




(8) 


, and dA = dxdy, A corresponds to the area in the (x,y) plane contained within 

^ A' r 

the region bounded by the closed curve ^,0 is the unit normal to the curve C 

<• C 'f 

in the (x,y) plane, is the coordinate velocity in direction X { 1. = \ ^ ) 

y ” - 

= Vj ) • Equations (6) 
dt 

to (8) represent the conservation of mass theorem in terms of area integrals 
in the (x,y) plane and line integrals evaluated on a curve c in the (x,y) 


and ^ is the coordinate velocity 
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plane. 


The curvilinecir effects are accounted for by the metrics h ,h , 

X y 

h and their derivatives. The term on the right-hand-side of Equation 
z 

(6) is ein elliptic source term and the second and third terms on the 
left-hand-side of Equation (8) are parabolic terms. 


6.0 


TRANSITION AI^D TURBUT RNCI-: 


In this section various models of tramsition cind turbulence are 
investigated and the proper ones are selected for incorporation into impeller 
computer code. Subsection 6.1 deals with the turbulent models, Subsection 
6.2 concerns the mixing length theory, emd separation is discussed in 
Subsection 6.3, and Subsection 6.4 considers transition. 

6.1 Turbulent Models 


With the advance of high-speed computers, turbulent flow problems 
have become amenable to numerical studies in the past decade. The develop- 
ment of turbulent models has contributed substantially to these studies. 
Though their progress is still in a preliminary stage, there is no shortage 
in the supply of models. The difficulty, from a user's point of view, is to 
select an appropriate model for his particular problem. All models of 
turbulence are supposed to be general and few cross comparisons between 
models are available. However, at the present time there is no definitive 
verdict as to the best turbulence model to employ. Thus a good rule in 
selection seems to be "the simpler the better". 

The adoption of models for turbulence naturally rules out the 
relatively more fundamental approach via statistical theory, which might be 
academically pleasing but unrealistic in engineering applications. In 
general, turbulence modelling is divided into two classes; those described 
by one algebric relation, such as the mixing length hypothesis, and those 
described by one or more differential equations governing some quantity 
like turbulence energy, turbulence vorticity or shearing stress. There 
are numerous examples in the latter class, generally referred to as the 
transport model, for example the classical Kolmogorov model (1942)^^ and 
recently the Saffman model (1970)^^ In adopting such a model, one must 
solve, in addition to the conservation laws, several differential equations 
from which turbulence stresses are determined. Limited success can be 
claimed in application of the tramsport models; they all seen to do quite 
Well in simple problems like turbulent boundary layers with small pressure 
gradients. 


; 
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Let us present herein the Saffman transport model for illustration. 
The model contains two variables: the energy density e and a pseudo-vorticity 

O » which cure assumed to satisfy the following non-linear diffusion equations. 





4 




(9) 




-y3/n,/ 


t 

Sj. 



( 10 ) 


where: t = time, 

Xj = Cartesian coordinates (j = 1,2,3) 

= mean density 

u^ = mean velocity components in the j-th direction 

yti = molecular viscosity coefficient 

S . . = mean rate of strain tensor 
13 


The numbers 0( , o(^ ' 


to be universal constants. 
= cr* = 1/2 
oc* = 0.3 
/i* = o/, 

5/3 ^ ^ 2 

t y _ 




^ = 


p* 


4 cr 


j * 

c/. 


-) 


2.5, based on experimental data. 


and Kj the Karman constant. 


<? 

c 


are assumed by the model 
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Equations (9) and (10) are integrated with an appropriate set of boundary 
conditions (which are not trivial) to yield e and i) . The eddy viscosity 
6 is related to e and by 


L- 


c 

n 


( 11 ) 


Saffinan’s model is but one of the many available schemes governed 
by two equations; some of the others are Chou (1945)^^, Harlow-Nakayama 
(1968)}'^ Jones-Launder (1972) Ng-Spalding (1972) etc. They all have a 
set of empirical constants, some even parametric functions. The complexity 
of the mathematical system and the uncertainty in those constants are 
inherent with all the models. Moreover, a set of non-linear diffusion 
equations generally introduces a new time scale in the computation, which 
is often substantially smaller than the convective or diffusive time scale 
for laminar type computation. The two-point boundary value problem also 
poses a tedious numerical task. However, the advcintage in this kind of 
turbulence modelling is also clear; they all attempt to depict the physics 
of turbulence transport, generation, dissipation and diffusion. In addition, 
some models (such as Saffman's) show the correct analytical behaviour near 
the wall (as demanded by the law of wall) . The predictive capabilities 
for incompressible boundary layer flows by those models are thoroughly 
established. Turbulent flows in more than two spatial dimensions, including 
separation, compressibility, rotational effects, and containing boundary 
layers interaction with shock waves have not been subject to examination 
by those models*. In short, the turbulence models, as promising as they 
are, have yet to be thoroughly tested by problems more complex than plane 
boundary layer flows. 

In view of the three dimensionality of the compressor problem, 
the desired economy in computation, and the added degree of complication in 

1 *7 

*Wilcox , applying Saffman's model, has shown good results in the study of 
turbulent boundary layer separation and reattachment at moderate (2.9C) 

Mach number . 
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the nonlinear equations, we must seek an alternative to the formulation by 
turbulence model equations. The alternative should be able to render a 
reasonably good description of the turbulent boundary layer development 
without a disproportional amount of computation time. 


6.2 


Mixing Length Theory 


The mixing length theory herein is the one originated by Prandtl 
amd subsequently modified by Van Driest^^, Cebeci^O other researchers 

to include the effect of compressibility. The formulation, in comparison 
with turbulence models, is quite simple. The Reynolds stress tensor** 

X is expressed by the eddy viscosity E . 


18 


t. 




(12) 


where ^ e = “ ' 2 ' 

The eddy viscosity E is then estimated by the mixing length theory which 
subdivides the sheaur layer into am inner and an outer region. Near the 
wall, we have 




by 


P 


(13) 


where; = 0*4 

y = normal direction from the solid wall 
u = velocity component parallel to the wall 
D - 1 - exp (^) 


A = 


Y 


**Total stress tensor 'T . . is given by . . = T. . 

ID ID ID 

the mean laminar stress tensor. 


with T. . being 
^iD ID 


-ly- 


I 


I 


= kinematic viscosity coefficient 
'lut - shearing stress at the wall 
‘^P/dx = pressure gradient in the direction parallel to the wall. 

In the outer region, the so-called Clauser defect law is used. 

( 14 ) 


where 


max 
* 

K 


k 


0.0168 

maximum value of u in the direction normal to the wall 
incompressible displacement thickness 



Jo 


c r • 

The upper limit d in the integral has to be defined to suit the compressor 

problem. It can be taken either as the mid-point of the channel defined by 

the blade-to-blade surface or the location where u occurs^. The eddy 

max 

viscosity £ is then evaluated by 


'f 


(15) 


A typical variation of (£ is sketched in Figure 5. When the flow 
field extends to infinity in the y-direction, £ is often multiplied by 
the Klebcuioff intermittency factor o , 



h 5.5 C-i")" _ 


( 16 ) 


which effectively makes decay away from the wall as it should by physical 

reasoning. The introduction of Y is not necessary in the compressor problem. 

The mixing length theory has enjoyed a large number of followers 

and many successes in applications. Cebeci eind Smitl?^ have api)li* d it 

successfully to incompressible cind compressible boundary layers, with 

and without separation-^ , with ma~s and heat transfei^O, as well as low 

Reynolds number turbulent flows^^*^^. Figure 6 shows a comparisori of 

22 

results obtained by the mixing length theory with measurements for a wasted 

+The definition of does not affect the value of li for a monatomic 
velocity profile. 


- 20 - 


T 


T 







.■'MODIFIED BY / 



0 T 


FIGURE 5, TYPICAL VARIATION OF EDDY VISCOSITY COEFFICIENT £ . 



FIGURE 6 : COMPARISON OF RESULTS OBTAILTD THROUGH THE MIXING LE!,'GTH TilEOHY 

AND MEASURI:MI.NTS of winter, SMITH AND KOTTA FOR A WAISITID BODY CF Kj;VOI.UTI UN . 
THE EXTERNAL KAQJ NUMBER M^ IS GIVEN, THE SKIN I WCTION COEFFl Cl EIIT C AND 
THE MOMENTUM 'IHICKNESS 6 A^ C0MPUT1;D A-ONG THE A>:TS. SOLID CUR'.T: Mi pj .i :."i . 
CEMECI'S FiiSULT WITH COKM'.CTION FOR THE TK/iNCVi:REE CURVATURE FMT'CT, DCTTED 
CURVE SIGNIFIES RIISULTS WITHOUT THE CORRECTION. THE EFFECT OF TrANCVUv.E 


CURVATURE IS NEGLIGIBLE IN THE PRJ:DICTI0N OF C 


M„o , AND THE REYNOLDS NUMBER R^, 
ALSO GI\T:N in THE FIGURE. 


f*. 


THE Fi%i:,r.STI<i AM M/TCH tU-'MiiEl- 


BASED ON THE UINGTH OF THE BODY L, AKl. 
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body of revolution. The calculated skin friction coefficient and 
momentum thickness 0 seem to agree with experiment quite well. Most 
recently, Cebcci* successfully applied tlie mixing length theory to tlio 
study of internal flows in area-dianging channels. However, the failure 
of the mixing length theory in the interaction problem of a strong shock 
with a hypersonic (free stream Mach nu ,er = 8.5) turbulent boundary layer 

2 D 

has also been reported by Baldwin cind MacCormack . It yields an incorrect 
pressure rise estimation on skin friction and heat transfer. However, 
the Saffmein predictions were also found to be inaccurate in the same problem 
by the same authors. One may conclude that the mixing length thnory , which 
is applicable to the attached boundary layer, is perhaps inadequate for 
the prediction of separated and reattached flows of the type examined by 
Baldwin and MacCormack. Again, tlie same conclusion is drawn for the 
transport model, Saffman's in this case. Fortunately, shock -boundary - 
layer interactions of this magnitude do not occur in the compressor 
problem. The track record of mixing length theory seems to certify its 
usefulness . 

In almost any discussion concerning the mixing length tiieory 
(or equivalently, the concept of eddy viscosity), the criticism that eddy 
viscosity should not be a local property inevitcibly arises. It is indeed 
true that turbulv. :e is a macroscopic phenomenon marked by eddies of 
finite size, and possesses a relaxation time similar to that r.f a visco- 
elastic solid. However, evidence has been gathered over the years to support 
the concept of a local eddy viscosity coefficient. It is believed th.it 
the mean-velocity gradient anu tlie turbulent shearing stress generally go 
up and down together, in particular, they go to zero together. The concept 
of eddy viscosity leads to accurate predictions of velocity profiles c'^cr 
a vast range of parametric inpiuts. 

2 7 

In summary, tlie advantages of the mixing length theory are that 

1. It is simple, requiring no addition.il differential equ.it i(<n 
to solve. Ihis is crucial in tlie comi'ressor calculation, since the existing 
routine is complex enough because of the geometry of the compressor. 
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2. It allows a realistic prediction to be made of the velocity 

2 8 

and the shear-stresses, and the general behaviour of boundary layer flows. 

3. Much experience in the use of mixing length theory has been 
accummulated and is available in the literature. 

Besides the comment on "local" property, the arguments agai..st 
mixing length theory are as follows: 

1. There is no successful evidence in predicting recirculating 
flows. However, the same comment applies to existing transport models. 

2. It implies that the effective viscosity vanishes where the 

velocity gradient is zero, ^ . Transport models, on the 

other hamd, do not provide a definitive relation between 6 and | j , 

they provide a set of differential equations whose solution presumbaly 
defines chat missing relation. 

3. The nixing length theory takes no account of the convection 
or diffusion of turbulence. 

The development of turbulence models is still in a preliminary 
stage; much modification to existing models is expected in the years ahead. 
In the absence of a clear-cut all purpose model, the one which has been 
experimented with the most, has shovv-n the most success, and which is 
simplest to use should receive first attention in our compressor studies. 
Hence, we selected the mixing length theory as our tool in turbulence 
studies . 


6.3 Separation 

In general, there are two types of boundary luyer separations: 
Iciminar and turbulent. Separations in compressors or diffusers may occur 
either way depending on the inlet conditions. In computational fluid 
dynamics, various criteria have been examined to identify the spearation 
point. 
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Separation on a plane or eoci symmetric flow is defined by the 

point where Uie wall shear vanishes, namely, 7^= 0. In lami nar flow, 

one can monitor the variation of tlio wall shear to locate separation. In 

addition, there are other simple criterion such as that based on the 

momentum integral metl»od of Tliwaites^^, and that of stratford^^ 'ough 

which laminar scjiaration is d-^fined. For cxiunple, laminar separate i is 

predicted when ) reaches a value of 0.102, where 

C is the local pressure coefficient and x the streamwiso distance from 
P 

the leading edge. We shall simply monitor the wall shear and the pressure 
distribution to pinpoint the laminar separation point. 

The prediction of turbulent separation is a much more difficult 
task. The current prediction methods can be divided into two grouj^s. Tiicse 
methods arc either of differential type (meaning that partial differential 
equations are solved) or of integral type (meaning that momentum integral 
or energy integral equations arc solved)^. The simplest integral method 
involves tiie monitoring of the shape factor H, H = S. /Q . ( S. is tl;e 

ill 

incompressible displacement thickness, B. the incompressible ir, omentum 
thickness). Wnen H reaches a certain value (H = 2.8 is used by .XcNally , 
and the range between 1.8 and 2.4 is frequently cited), then se,)aration of 
the turbulent flow is assumed to occur. Again, we shall just monitor t)ie 
wall shear and extrapolate it to zero for the location of the separation 
point. In past computations, the peak of the wall pressure and the vanishing 
of the wall shear locate separation jointly. 

In short, we -.hall introduce no additional criterion for the 
determination of laminar or turbulent separation other than its natural 
definition through the vanishing of wall sliear. Since we solve the full 
Navic r-Stokes equatjons numerically, we do not )iave to change gover.ning 
equations after sej'aration takes place. 


+T)ic 1008 Al'OSR-lI-T-Standard Conference on "Com]>utat ic<n of Itinnul e.nt Houi.d.iry 
Layers" provides a critical cv.jluation of tJiose mothods^^. 
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RADIAL IMPELLER PROBLEM 

To fully describe the debug problem selected, which was designated 
as "Problem 1.0", nine items are discussed. They are (1) input conditions, 
(2) geometry in cylindrical coordinates, (3) curvilinear coordinate 
system, (4) boundary conditions, (5) meshes, ( 6 ) inviscid solution, (7) 
medium viscosity, ( 8 ) constant speed, , at which the blade-to-blade 

surfaces of calculation move from the hub-to-the-shroud, and finally, the 
results (9) of the numerical computation for problem 1.0. 

7.1 Input Conditions 

The radial impeller which served as the debug problem was selected 
by Dr. T. Katsansis of NASA Lewis Research Center. Input flow properties 
for the problem are as follows; 

laboratory inlet total temperature = 536 °R 
laboratory inlet total pressure = 861 psfa 
rotational speed = 4031.70 rad/sec (38,600 rpm) 
specific heat ratio = 1.667 
gas constant = 38.73 ft/Oj^ 

The above specific heat ratio and gas constant cure for Argon. 

7.2 Gecanetry 

The radial impeller geometry is presented in the cylindrical 
coordinates r, 0, amd X 3 . These coordinates are defined in Figure 7 in 
terms of Cartesian cooru '.nates Xi, X 2 , and X 3 . The sense of the angular 
rotation is also indicated. Ic is seen from Figure 7 that a left-handed 
coordinate system is being employed. Figure 8 presents traces of the hub 
and shroud lines in a half-plane through the axis of the impeller, i.e., a 
meridional plane. The hub and shroud radial coordinates are presented as 
functions of axial distance. Solid lines indicate the accual geometry of 
tne machine, while dashed lines indicate formula approximations. The hub 
is approximated by an ellipse with its major axis on the radial axis, and 
the shroud is approximated by a super-ellipse having its major axis in the 
axial direction. The elliptic and super-elliptic formulas are also shown 
Figure 8 . dashed hub shroud liner extend above the discharge, 

so that the region of calculation contains the discharge. The hub elliptic 
formula approximation produces an inlet area about 9% greacer than th. 
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actual inlet area. Although approximate, the hub formula simplified 
development of the IFFC computer code. The angular coordinates of the 
pressure and suction blade surfaces are shown in Figure 9 as functions of 
eixial distance X^. The solid lines indicate the traces of the blades on 
the hub, while the dashed lines indicate shroud blade tiaces. 

7.3 Curvilinear Coordinate System 

An eucisymmetric orthogonal, curvilinear coordinate system was 
used to solve debug Problem 1.0. Consider the curvilinear coordinates 
x,y, and z. The surfaces x = constant were selected as half-planes through 
the axis of rotation of the machine, i.e., meridional planes. The surfaces 
y = consteint and z = consteuit were obtained by rotating two orthogonal curves 
in the meridional plane about the axis of rotation of the machine. Since 
the hub was approximated as an ellipse, elliptic coordinates were used to 
establish the y and z surfaces. A family of confocal ellipses defined the 
z surfaces and a family of hyperbolas, orthogonal to the ellipses, defined 
the y-surfaces. Consider the interior ellipse shown in Figure 8 . Tliis 
ellipse is labelled with the constant value z =-1.196. The value of z 
is determined from the following formula: 

tanh z =“ 7 ( 2 ) 

where B and A are the lengths of the minor and major axis, respectively, 
of the interior ellipse. The negative sign permits z to increase as the 
blade-to-blade surface moves from the hub to the shroud. The orthogonal 
hyperbola to this ellipse is labelled witii the constant value y. The value 
y is defined as the single the asymptote of the hyperbola makes with the 
radial axis. 

The transformation equations from cylindrical r, 6 , X 3 space to 
curvilinear x,y, z space are shown below; 

r = Rq= C Cosh Z Cos y 


X 3 = C Sinh Z Sin y 

where Rq is the maxium radius of the impeller and C is the focus of the 
elliptic and hyperbolic coordinates. Formulas for the metrics of trans- 




formation Equations (3),i.e., hx,hy,hz are derived in Table I as well as 
their derivatives. Application of Equation (3) to the geometry of Figures 
8 and 8 results in the transformed geometry in the x,y,z space. The trans- 
formed geometry is shown in Figures 10 and 11. Figure 10 presents the hub and 
shroud lines in the y,z, plane. The hub is a horizontal straight line, since it 
is an ellipse. The shroud is still a curve, since it is a super-ellipse 
with a reversal of major and minor axis with respect to the hub ellipse. 

The interior ellipse of Figure 8, corresponding to z = -1.196 radians, is 
shown as a dashed horizontal straight line. Figure 11 shows the pressure 
euid suction blade surfaces in the (x,y) plane. The solid lines indicate 
the blade traces on the hub, while the dashed lines indicate the blade 
traces on the shroud. In curvilinear space the calculation will take 
place in (x,y) planes which move from the hub to the shroud as the 
parameter z increases. 

7.4 Bo\xndary Conditions 

Boundary conditions for the impeller problem in the (x,y ) planes 
of calculation are indicated in Figure 12. At the upstream boundary 
of the region of calculation uniform inlet conditions are specified. 

Along the pressure and suction blade surfaces no slip flow is enforced. 

At the lateral boundaries upstream of the inducer and downstream of the 
discharge periodic boundary conditions are enforced. Finally, the back 
pressure is specified at the downstream boundary of the region of calcu- 
lation. 

To expedite development of the IFFC computer code, the upstream 
boundary was placed at the inducer in the solution of debug Problem 1.0. 

Inviscid conditions, discussed in Section 7.6, were prescribed along the 
upstream boundary. However, it is emphasized that boundary conditions of 
Figure 12 will be employed in the solution of all problems subsequent to 
the debug problem. 

7 . 5 Meshes 

From Figure 11 it is seen that the (x,y) blade-to-blade surface 
distorts as z increases from its value at the hub of z —1.22524 radians 
to z= -1.14 radians near the shroud. Since the blade shape in the (x,y) 
plane distorts with z, the finite difference mesh must distort as well . 

Thus, a subroutine was developed to automatically distort the iinite 
difference mesh in accordance with the blade geometry. Two meshes, 
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TABLE I: VALUES /-ND DERIVATIVES OF THE TRANSFORMATION 

SYSTEM OF PROBLEM 1.0. 
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corresponding to z ■ - 1.22524 radians and z •■-1.14 radians are shown in 
Figvires 13 and 14| respectively. Each mesh is formed by the intersection 
of 20 streamline-like-lines and 39 potential**like-lines > i.e.i 780 points. 
The streamline-like lines are spaced closer in the vicinity of the blades 
than in the center of the blade passage. 


7.6 


Inviscid Solution 


In order to solve the equations of motion shovm in Appendix A, 
the inviscid solution is required. The inviscid flow field serves two 
purposes. First, the inviscid field at the hub provides initial conditions 
for the viscous calculation: for debug Problem 1.0 the hub boundary layer 
was neglected. Second, as discussed in Section 4.0 the inviscid flow 
field is the zeroth iterate in the interation procedure. 

The inviscid flow field for debug Problem 1.0 was solved for by 
Vanco (Ref. 36) using the meridional velocity gradient method of Katscuiis 
(Ref. 10) . The velocity vector eind pressure fields were calculated on the mean 
mean hub- to- shroud stream surface of the impeller channel. These proper- 
ties were specified along five streeunlines, as well as the suction and pressure 
blade velocities associated with each r, X3 point along the streamlines. 

At a given mesh point in the flow field the velocity vector, 
specific internal energy, pressure, £ind density were determined in the 
following manner. First the velocity vector was found by linear inter- 
polation in the inviscid field of Vanco. The rothalpy, H = E +w2/2-r^u<»^/2 
which is invariant along inviscid streaunlines , was then used to compute the 
specific internal energy in terms of the velocity and radius at the given 
point, i.e.. 


1 - 


w2 CL»2y2 
2 


(4) 


where E^is the inlet stagnation specific internal energy in the laboratory 
frame, W is the magnitude of the relative velocity vector, r is the local 
radius, cu is the amgular velocity and /is the specific heat ratio. 

Pressure was calculated from the given meam stream surface pressure by 
assuming insetropic flow along each blade-to-blade circular arc associated 
with r euid X^. The density was then determined from the equation of state. 
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The flow field above the discharge was computed from relations 
governing flow in a vaneless diffis'jr. The mass# angular momentum in a 
laboratory frame* and rothalpy were conserved. The viscous mixing above 
the discharge was partially accoiinved for through an entropy gain; the 
density at the downstream boundary was reduced to 95% of it's isentropic 
value. 

The inviscid relative velocity field at the hub* z • - 1.22524 
radians* is shown in Fig. 15. These vectors have magnitudes proportional 
to the particle velocities in the (x,y) plane; their tails emanate from 
the mesh points of Figure 13 . The e symbols indicate the pressure and 

t 

suction blade surfaces. A value of ^/z was added to the angular x values 
of Figure 13 to produce positive ordinate values. The locations of the 
pressure and suction blade surfaces are reversed between Figures 13 ^Lnd 15 
because the abscissa of Figure 15 is located on the top of the page* while 
the abscissa of Figure 13 is located on the bottom of the page. It is 
seen from Figure 15 that the velocity profile is linear between the 
pressure and suction blades. Furthermore* Mach number calculations in 
the vicinity of the inducer indicate transonic flow. For example* at the 
upstream boundary of the region of calculation* i.e.* the inducer* the 
suction blade Mach number is .95. 

7.7 Medium Viscosity 

The meshes of Figures 13 and 14 have zone widths in the neighbor- 
hood of the suction and pressure blades which are too coarse to define 
a thin boundary layer. For example consider the inviscid field at the hub 
shown in Figure 15. At the inducer the suction blade Reynolds number per 
foot for Argon is 2.09 x 10^. The hub ellipse is .26 feet in arc length* 
so based on this dimension the exit Reynolds number for Argon is about 
545*000. This Reynolds numbers would probably produce a turbulent boundary layer 
thinner than the width of the first layer of zones adjacent to the suction 
blade surface. Therefore* to resolve the boundary layer with the meshes 
of Figures 13 and 14* a Reynolds number reduction is required. 

Flat plate analysis indicated that for an inducer suction blade 
Reynolds number of 20*000 per foot* the leuninar boundary layer at discharge 
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is contained in about five zones of the mesh. Based on the arc length 
of .26 feet of the h\ib ellipse, the discharge Reynolds number is then 
5000. 

Therefore, in order to get meaningful results with the meshes 
of Figures 13 and 14 it was necessary to invent a fictious medium having 
the compressibility properties of Argon cuid the viscous properties 
appropriate to a Reynolds number of 5000. Thus, the flow field of 
debug problem 1.0 is quite far from actual iiq>eller flows and can only be 
considered in a qualitative sense. 

7.8 Specification of Uqq i Speed at which Blade to Blade Surfaces 

of Calculation Move from Hub-to-Shroud 

To run debug Problem 1.0, the speed Ucs at which the (x,y) 
blade-to-blade planes move from the hub to the shroud must be specified. 

The final steady solution will be independent of Uoo » however, intermediate 
solutions will depend on the magnitude of Uc© • The speed Uqq must be 
small enough to permit viscous diffusion effects at the blade surfaces to 
build up boundary layers which subsequently separate. The time it takes 
a particle at the inducer to travel to the discharge is the characteristic 
time, t , for boundary layer build-up (Reference 3) . 

The speed U^q was determined in terms of the characteristic 
time t. in the following manner, consider flow along the suction blade at 
the hub. The average velocity is adiout 583 fps and the hub arc length is 
.26 feet. Therefore, the characteristic time is .444 ms. Since the 
time it takes the boundary layer to develop is t, let us assume that another 
characteristic time is necessary to permit the boundary layer to separate. 
Hence, approximately 2t characteristic times are required for the (x,y) 
blade-to-blade plcuie to go from the hub to the shroud. If 1.9 character- 
istic time passes in the time period that the (x,y) blade-to-blade surface 
moves from the hub to the shroud, then the appropriate speed is U^ = 

21.8 fps. 



7.9 


Radial Impeller Numerical Results 


I Prc^lem 1.0 was run 750 cycles*, or Icmg enough in time 

i 

t for the z parameter to increase frc»n z * -1.22524 radians at the hub to 

j 

i z -1.194 radians; the (x,y) blade-to-blade plane moved approximately 37% 

of the total increment in z between the hub and the shroud at discharge. 

I This partial solution of Problem 1.0 required 23 minutes on the CDC 7600 

coiq>uter. 

Shortly after the calculation commenced, it was foiuid that there 
was a mass imbalance in the initial conditions. The elliptic formula 
approximation to the hub geometry (see Section 7.2) increased the inducer 
flow area by 9% from the actual flow area. Vanco's inviscid solution was 
not corrected for this geometry change. Therefore, a mass imbalance was 
produced in the initial conditions of Figure 15; more mass flux entered the 
system than exited from the system. 

The okass imbalance is indicated in pressure distributions along 
the blade surfaces. In Figure 16 pressure distributions on the pressure 
blade surface are shown for three values of the z parameter. The 
abscissa of Figure 16 represents the eingular coordinate y, while the 
ordinate is the ratio of the local to inlet stagnation pressure, i.e., 
p/p^ . Curve 1 represents the initial pressure blade distribution at 
z » -1.22524 radians; the initial exit pressure ratio is 1.635. Curve 2 
represents the distribution 100 cycles after the start of calculation, 
i.e., z = -1.2189 radicUis. There is a pressure peak in the center of the 
chamnel followed by a deep reire faction in the radial portion of the 
in^ller. The deep rarefaction and fixed high exit pressure induced back 
flow at the downstreaun boundary, a condition of impeller surge. To prevent 
surge the back pressure was lowered to a ratio of 1.47. At the lower 
back pressure level outflow was maintained at tl.e downstream boundary and 
the problem was continued. 

*A cycle of calculation consists of updating the dependent variables 
of motion through one timestep over all the mesh points. 
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The pressure wave and trailing rarefaction calculated at cycle 
100 started moving towards the inducer and the solution in the radial 
portion of the impeller converged. Convergence was demonstrated by 
inspection of the flow field over a small change in the z coordinate. For 
a given small change ir z the flew field in the radial portion of the impeller 
ch6uiged slightly, while flow near the inducer changed markedly. Curve 
3 of Figure 16, which corresponds to cycle 550 or z =-1.196 radieins, is 
converged in the rad: al portion of the in^eller. The solid line represents 
the converged region of the flow and the dashed line represents the un- 
converged region. 

The upstream moving pressure wave finally impacts the upstream 
boundary, reflects from it, and amplifies. The inviscid conditions 
prescribed at the upst; earn boundary cause reflection and aunplification 
of the pressure wave. This saune phenomenon was observed in previous 
cylinder calculations started from impulsive initial conditions (Ref .2} . 
Problem 1.0 cor'd not be continued beyond this point without moving the 
upstream boundairy upstream of the inducer. 

The three-dimensional flow field in the radial portion of the 
impeller has converged amd is very interesting. Results aure presented 
for the radial portion of the imp'^ller in the remainder of this section. 

The sequence of events as the flow develops in the impeller 
chamnel is illustrated in the velocity vector plots of Figures 17 and 
18. Figure 17 shows a velocity vector plot at cycle 100 (z = -1.2189 
radiauis) amd Figure 18 shows the velocity field at cycle 450 (z = -1.991 
radians) . In Figure 17 bovndary layers are seen on both the pressure amd 
suction blade surfaces. A flow instaUsility is beginning to occur at the 
downstream end of the suction blade surface. At cycle 450 (z = -1.991 
radiains) the flow has converged in the radial portion of the impeller. 

A thicker pressure blade surface boundary layer thain present at cycle 100 
is clearly indicated in Figure 18. Furthermore, the suction blade boundary 
layer has separated amd a large vortex is in evidence on the suction bla'^e 
surface near discharge. The vortex takes up almost half the channel 
width between blades. The reduced chamnel flow area causes an acceleration 
of the flow about the vortex; large vectors are in evidence just above the 
vortex. This large vortex is consistent with the size of vortices previously 
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FIGURE 18- RELATIVE VELOCITY VECTOR PLOT OF THE VISCOUS IMPEUXR PUR* FIELD OR TOE ELLIP1 
BIAOC-TO-BIAOE SURFACE OF REVOLUTION AT a » -1.1991 RADlANSi ♦ SYrCOLS INDICA': 
TOE PRESSURE AND SUCTION BLADE SURFACES } StmFAO: FOVED 32* FROM HUB ATdISQUL<I 
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determined cibout cylinders at lew Reynolds number (Refer-inces 1 and 2) . 

A cosqparison of viscous and inviscid velocity fields is pre- 
sented in Figures 19 and 20 on an (x,y) plane of calculation which has 
moved about 33% of the total increment in z at the discharge. In the 
viscous flow field of Figure 19, which corresponds to cycle 550 
(z = -1.196 radians), the separated region has grown larger and feeds 
into the boundaury layer on the suction surface. The subsonic nature 
of the flow causes the suction blade velocity vectors upstream of the 
separation to adjust to the vortex. Figure 20 shows the corresponding 
inviscid flow field at z = -1.1960 radians. Due to the absence of viscosity, 
the inviscid suction blade flow does not separate. 

The velocity field in the radial portion of the impeller at 
the final cycle calculated, i.e., cycle 750, is shown in Figure 21. A 
well formed vortex is seen in Figure 21 which extends aft of the discharge. 

A reversed flow profile is clearly seen in Figure 21 . The results of 
Figure 21 indicate that the velocity field is highly non-uniform above 
the discharge plane. 

A comparison of the viscous and inviscid pressure distributions 
in the radial portion of the impeller is presented in Figure 22. These 
data correspond to z = - 1.196 radians or on an (x,y) blade- to-blade 
plane of calculation which has moved 33% of the total z increment between 
the hub and shroud at discharge. The pressure surface comparison (Figure 
22b) indicates that the viscous pressures are no more than 8% high them i 
the inviscid pressures upstream of the station in the channel where the 
bac)c-pressure influences the discharge flow. The lower bac)c -pressure in 
the viscous suction blade rurface pressure of Figure 22a drops off at 
the separation point to nearly coincide with the inviscid solution. The 
rapid drop in pressure in the viscous case is consistent with the in- 
crease in the flow velocity just cUx>ve the vortex (see figure 21) . 

Although the flow field in the inducer region has not yet 
converged, it is clear that the IFFC computer code has duplicated, 
at least qualitatively, the flow phenomena that have been observed to occur 
in the radial portion of an impeller (Ref. 37). We therefore conclude that 
the IFFC computer code works for laminar flows. 
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FIGURE 22b COMPARISON OF VISCOUS An15' PRESSURE VARIATIONS ON THE 
RADIAL PORTION OF TOE IMPELLER; z = -1.1960; PRESSURE SURFACE. 
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FIGURE 22a COMPARISON OF VISCOUS AND INVISCID PRESSURE VARIATIONS ON THE 
RADIAL PORTION OF THE IMI’ELLER; z = -1.1960; SUCTION SUWACT:. 



8.0 BACKSWEPT IMPELLER PROBUIM 

i 

A proof of principle impeller problem was run to check out the 
turbulence model which had been incorporated in the IFFC code. The impeller 
geometry^ which was selected by Dr. T. Katsanis of the NASA Lewis Research 
Center, was that of an advanced backswept compressor developed by CREARC, Inc. 
The basic geometry and design operating conditions of the backswept impeller 
problem were as follows: 

Rotational speed , 7S0Q0 PPM 

Tip diameter 15.95 Cm. (6.26 in) 

Design pressure ratio 8:1 

2 

Inlet total pressure 2117 Ib/ft 

I Inlet total temperature 519 ® R 

Impeller tip speed 2055 ft/sec 

Discharge Reynolds Number* 1.43 x 10® 

The rotor geometry is presented in the cylindrical coordinates 
T, 6 , Xj in Figures 23 and 24. Figure 23 presents traces of the hub and 
shroud lines in a half-plane through the axis of the impeller, i.e., a 
meridional pl^me. The hub and shroud radial coordinates are presented as a 
function of axial distance. The hub and shroud lines extend upstream of 
the inducer and eibove the discharge so that the region of calculation 
contains them both. The angular coordinates of the pressure and suction 
blade surfaces are shown in Figure 24 as functions of eucial distance X^. 

The solid lines indicate traces of the blades on the hub, while the dashed 
lines indicate shroud blade traces. The regions upstream of the inducer 
and downstream of the hub are also indicated in Figure 24. 

A fine finite difference mesh was incorporated in order to 
f adequately define the boundary layer. The mesh consisted of 30 J-lines 

(streamline-line) emd 101 K-lines (potential-like) . A grating factor 
(g) of 1 C6 was usod to space the J-lines. The hub plane mesli is illustrated 
I in Figure 25 and the blade-to-blade surface 23% of the distance between 

I hub and shroud is shown in Figure 26. 

The invicid flow solution for the backswept impeller was solved 
1 using the meridional finite difference method of Reference 39. The hub 

^ inviscid flow field, which serves as the zeroth iterate for the viscous 

, calculation, is shown in Figure 27. It is noted that at the inducer 

inlet the relative velocity is roughly the same on both the blade 
pressure and suction surfaces. The velocity profile remains relatively 

•Reynolds number based on an average inviscid relative velocity along the 
hub and distance along the hub. 
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constant until approxim.itely half way through the channel, when the suction 
blade velocity becomes significantly larger than the pressure surface velocity. 
This result is not in agreement with observed flow phenomena for centrifugal 
impellers. That is, the low-flow region is observed to occur at the suction 
blade surface near the discharge. 

Substantially different results were obtained for the viscous 
solution to the backswept impeller problem. Relative velocity plots for 
the blade-to-blade surface 19% of the distcuice from the hub to shroud are 
shown in Figures 28-30. The inducer region, which is shown in Figure 29, 
has a slight separation on the blade suction surface near the inlet. At 
the discharge the flow velocities near the suction and pressure surfaces 
and across the channel are nearly equal, whereas the inviscid calculation 
(Figure 27) predicted very low velocities on the pressure surface. There 
is no indication of a suction surface separation at the discharge like 
that obtained in the radial impeller problem. No conclusion can be drawn 
with regard to the effectiveness of backsweep in reducing or eliminating 
flow separation because the radial impeller was calculated for a very low 
Reynolds number (5000) with leoiinar flow, whereas the backswept impeller 
was calculated at a high Reynolds number (1.43 x 10^) with the turbulence 
model operational . 

The viscous solution for the surface 72% of the distance between 
hub euid shroud produced relative velocity profiles shown in Figures 31-33. 

The results were similar to the 19% surface except that the velocity 
gradient near the suction blade surface was not so pronounced. Also, 
except at the inlet, the velocity profiles were relatively uniform from 
blade-to-blade and continued to be so all the way to the discharge. 

The solution for the 98% surface, shown in Figures 34 and 35 
produced results which were similar to the 72% surface calculation. 

Likewise, the 99.95% surface, which is shown in Figure 36, looks much 
like the 98% surface, although the 99.95% surface is located in the tip 
cleareuice region. There is no blade to influence flow at this surface, 
but the relative velocity plot shows that the blade viscous effects are 
felt despite this fact. Indeed, the relative velocity (boundary layer) 
profile near the hypothetical blade surface is much like one would expect 
if a blade were present. 

The static pressure contour plot for the 98% blade-to-blade 
surface (z = .22451 radians) is shown in Figure 37. The static pressure 
ratio shown is referenced to the inlet stagnation pressure (P^) . There 
is a region of low static pressure near the inducer inlet (P/P^ = *8) on 
the suction side of the channel, which indicates a region of accelerated 
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flow due to the airfoil. The inducer or axial-flow portion of the channel 
decelerates the flow steadily until there is a static pressure ratio of 
about 1.4 as flow begins to enter the radial portion of the rotor. The 
rate of static pressure rise thereafter is increased rapidly because of 
the rotor centrifugal energy input. Finally, at the discharge an average 
chemnel static pressure ratio of about 4.6:1 is achieved cuid there is a 
relatively uniform profile across the channel. There is no separation 
indicated because the static pressure rise continues throughout the radial 
portion of the flowpath. This is attributed to the stabilizing influence 
of the backswept blading. 

Very similar results are illustrated for the 99.95% surface 
(Figure 38), which is located in the tip clearance region at the shroud. 

In Figure 39, hub, 19% and 77% surface calculations of suction surface 
static pressure ratio as a function of Y (axial coordinate) are 

shown. These results also indicate the absence of separated flow because 
diffusion continues all the way to the discharge. 

Relative velocity ratio contour plots are presented in Figures 
40-43 for the blade-to-blade surfaces located 19%, 72%, 92%, and 99.95% 
of the distance from the h\ib to the shroud. The relative velocity ratio 
V/V*, is the ratio of flow velocity to critical flow velocity according 
to the following*: t 


V/V* = M 


X-H 

2+()f-l)M^ 


At the 19% surface (Figure 40) the contour plot indicates roughly sonic 
relative velocities at the inducer inlet with diffusion down to a velocity 
ratio of .4 - .6 at the impeller discharge. The cont-» ir plot for the 92% 
surface, which is shown in Figure 42, suggests that most of the flow at 
the inducer tip is at a relative velocity ratio of 1.2 or higher. A 
smooth, even diffusion rate is indicated, and a relative velocity ratio 
of roughly .6 is obtained at the discharge. 

Contour plots of relative total pressure ratio are shown in 
Figures 44-46. Relative total pressure ratio is defined as the 
stagnation pressure calculated at a given point divided by the ideal 
stagnation pressure which would have occured if the process were isentropic. 

*Where M is defined as the ratio of the local relative velocity and the 
local sound speed. 
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The contour plot of the 19% surface, shown in Figure 44, indicates that for 
nvost of the flow the relative total pressure begins to deviate from ideal 
conditions at adaout 40% of the meridional distance through the passage. At 
the discharge a value of the relative totail pressure ratio of roughly .92 - 
.96 is calculated for flow in the mid-passage. The regions of low total 
pressure are observed to be near the suction and pressure blade surfaces. 
Similar results eure indicated for the 72% and 92% surfaces shown in Figures 
45 and 46. At the 98% surface it is observed that a region of relatively 
high (.96 - .98) total pressure ratio occurs at the discharge in the mid 
passage region. Again the main flow losses are calculated to occur near 
the blade boundaries. 

An independent boundary layer computation, using Mager's integral 
turbulent boundary layer analysis , was performed to verify the numerical 
separation characteristic. The calculated suction blade momentum thickness 
(§) in the boundary layer for the 19% and 77% surface calculations is shown 
in Figures 47 and 48. Since the ratio of (R) to the initial value (g) ^ stays 
relatively constant downstream of the initial inducer portion of the channel, 
there is no separation and no reduction of static pressure rise capability. 
Nowhere does the suction blade momuntum thickness go below its initial value; 
hence, turbulent boundary layer theory indicates that the suction blade flow 
should ;.ot separate. 

Distributions of the ratio of eddy viscosity to the molecular 
viscosity along the suction blade surface are presented in Figure 49. 

The boundary layer turbulence is generally increasing along the blade, 
especially in the radial flow region prior to discharge. 

The viscous calculation of the backswept impeller flow field with 
the IFFC blade-to-blade computer program ran 8073 cycles or surfaces from 
hvib to shroud. The 30 x 101 mesh consisted of 3030 zones and required 9.3 
hours of computer time on the CDC 7600 computer to complete the problem. 

It is recommended that additional efforts be made to reduce the comput- 
ational time requirement to make this solution a practical tool for 
utilization by the compressor aerodynamicist. 
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VISOOSI'n' TO MOLECULAR VISCCSITi’ ALONG THE SUCTION 


9.0 


CONCLUDING REMARKS 


A viscous blade-to-blade computer code for computing the flow 
field in centrifugal impellers has been sucessfully developed. The program 
was used to calculate the flow field of both a radial and backswept com- 
pressor impeller. Whereas the radial impeller problem indicated a large 
separation region on the suction blade surface near the discharge, the 
backswept impeller calculation, which included the mixing length turbulence 
model, did not separate. Mo conclusions can be drawn with regard to the 
effectiveness of backswept blading in reducing or eliminating flow separation, 
because the radial impeller was calculated for laminar flow at very low 
Reynolds number (5000) , and the backswept impeller was calculated at a high 
Reynolds number (1.43 x 10^) with the turbulence model included. 

Bie backswept impeller problem requires 9.3 hours of computer 
time on the CDC 7600. It is recommended that an effort be made to improve 
calculation efficiency to reduce the costs associated with the blade-to-blade 
solution. Also, it is recommended that the IFFC code be modified to provide 
the additional capability of calculating the flow in cross-section com- 
putational planes. This modification will enable the program to calculate 
shroud viscous effects, tip clearance effects, and comer vortices. 
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Appendix A 

DEVELOPMENT OF EQUATIONS OF MOTION IN ROTATING 
QRTHOGCWAL CURVILINEAR COORDINATES 


The coordinate system upon which the iteration t 2 d(es place controls 
convergence of the calculational procedure. Since flow is 

confined to an impeller blade channel, a generalized coordinate system, 
whose cixis follows the channel geometry, will be utilized to converge 
the iteration as rapidly as possible. Consider the generalized coordinates 
(x,y,z) shown in FigureAl; the surface x = constant is a mid-channel surface, 
surface y = constant is a blade-to-blade surface, and the surface z = constant 
is an orthogonal channel surface. The transformation of cartesian equations 
in coordinates, X^, X^, X^ to the generalized curvilinear coordinates x,y,z 
is presented in the following paragraphs. The development takes place in 
the following three steps: 

(1) The metrics and generalized basis vectors are derived. 

(2) Coriolis and centrifugal acceleration terms are developed 
in generalized coordinates. 

(3) The generalized equations of motion are presented. 

The rotating cartesian coordinates X^, X^, X^ are related to the 
orthogonal generalized coordinates x, y, z as follows: 


= fj^ (x,y,z) ^ 

^2 ^2 

X 3 = f 3 (x,y,z) ^ 


(Al) 


The metrics and basis vectors of the transformation can be determi.ned from 
Equations (Al) . An element of a length, ds, in cartesian coordinates is 
expressed in generalized coordinates as follows: 


.2 ^^2 ^ 2,2 ^ 2^2 
ds = h dx + b dy + h dz 

X y * 


(A2) 
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FIGUREiil: Orthogonal surfaces in the channel of a centrifugal ir.pcxler 

which define the curvilinear coordinates x, y, and z. 



where : h - 

X 


ft * 


The parameters, h^, h^, and h^ cure the metrics of the transformation. The 
vinit basis vectors of the generalized coordinates are related to the 


basis vectors , and 

as follows ; 

i = h 

— X 



j = h 
^ y 



k = h / 

^ 1, 4 3^ • ' -) 

ixx 

z ( 



The Coriolis and centrifugal accelerations, in the directions, x,y 
and z can be determined from a scalar product of their cartesian components 
with the basis vectors of Equations (A3) . The Coriolis accelerations in 
directions j_, and are, respectively 






where u, v, and w are the components of velocity vector u in the i_, j_, and 
Indirections, respectively, oj is the angular frequency, and CO is the angular 
frequency vector pointing in the direction. 
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The centrifugal accelerations in directions k» are, 

respectively 

Wxx;-i = (A7) 

('1.1) y tf/D-J = (A8) 

((^y Wv (A9) 

where r is the position vector in cartesian rotating coordinates X^, X^, and 
^ 3 - 

The Coriolis and centrifugal acceleration components of Equations 
(A4) to (A9) are added to the Eulerian set of equations of motion in 
generalized orthogonal coordinates (Ref . 36) • The final relations are as 
follows : 

Continuity 


( f wj =0 


(AlO) 


where ; 


(A - r 4 / j* 4- vv ^ 




and ^ is the density. 
X - Momentum 
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where ; 


y - Momentum 


^ i. ^ J + > 




where: 


z - Momentum 


^y - ^e? ~ '*' J 




<h(ii,t 'ifi? 5itA - 5i.«'i!!' . £?- /y-1. 


(A13) 


where : 


- 2?ii f J ^ C?'# y^> 

Internal l:nergy Equation 


^ ‘' ”1 


■h ^ 


^(K)J ^ )J 


(AM; 
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where E is the specific internal energy, j 

normal stress components, and shear stress 

com£>onents . 
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Appendix B 


CONTINUITY EQUATION IN GLUERALIZED COORDINATES 

In this appendix continuity Equation (4), in (x,y,t) space, 
is transformed to generalized coordinates The transformation 

equations between x, y, t space and the generalized coordinate s ,7 ,^) 
arc presented, and from these relations the integral equation for mass 
conservation is derived. 

Equation (4) is written in terms of the Eulerian coordinates 
(x,y,t). In the calculation the trace of the boundary of the impeller 
channel in the blade-to-blade surface must distort with time. 'I’lierefore, 
the continuity equation must be formulated in a generalized coordinate- 
system^, , X". The generalized coordinates arc defined as follows: 

t = T: (HI) 

X == f ( ^ , 7 ) (H2) 

y = g(^ ,1 ) (R3) 

and f ,7 .0) ,g()" ,7 ,0) = 7 , whore f- = > J - 1 , Ua - ) = i. 


Equation (Rl) is differentiated with re: pect to t, x, and y, 
respectively. Tin- result is as follows: 


O^J: 

o4r 



o'X 

o>X 



ll 


^ 0 J 


(U4) 


Di f 1 erontiati on of Equations (H2) and (R3) results in the followinu: 


■ hW 

Ml 

»x ■ [fjg^ - gjf^] 


»y Cfjg^ - gjf^] J 

Consider the function G(^,7 /T* ) • The derivatives of this function are as follows: 
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Equations (B3) - (B6) produced Equations (B7) , (E9; , and (B9) . 

Using Equations (B7) , (B8) , and (B9) , the continuity Equation (4) 

is transformed to the generalized coordinates ( r ■ >?.T ) . The trans- 


formed relation is as follows; 




1 


I 


where ( and ( define differentiation with respect to x and y, 
respectively, eind the grid velocity components and are defined in 
terms of deviatives of the curvilinear coordinates x and y. 


S = 


efZ 


dr 


= f. 


d-1. = 


(BID 


The symbol J represents the Jacobian of the trams formation, i.e,, 

1 


(X, y, t ) 




where ; 


dxdy = J^d^ d*7 


(B12) 


(B13) 


Equation (BIO) is multiplied by the area increment d^ d^ amd the resultant 
relation is the final continuity equation as follov;s; 




where dA = dx dy, A corresponds to the area in the x,y plane contained 

/s 

within the region bounded by the closed curve C, n is the unit normal to 
the curve C, £ is the particle velocity vector in the (x,y) plane as 
defined by Equation (7) , and is the coordinate velocity vector in tl.e 
(x,y) plane as defined by Equation (8). In the integration process 
use was made of Equation (B13) to convert integrals in d^ 0*1 to i r.^-cgral:; i:i 
dxdy. Furthermore, Leibniz's rule was used to permute differentiation and 
integration and Gauss's theorem was used to convert area integrals to linc' 
integrals in the (x,y) plane. 



